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(^^BSTRACT: numerical simulations are used to investigate the origins of the different wave velocities measured 
04n dense granular samples assembled with different methods. Glass bead packings are prepared in the lab either 

I i^y pouring and vibrating the dry material in a container, or by mixing with a very small amount of a viscous 

lubricant. Lubricated samples, although less dense, exhibit significantly higher wave velocities for confining 
C^-essures in the 100 kPa range. Numerical predictions for elastic moduli agree much better with experimental 
c/Jesults when the computational preparation of the samples mimics the laboratory one, albeit in a simplified 
^^anner. A plausible explanation to the laboratory observations is that the coordination number, which influences 
'^he material stiffness more than its density, is notably higher in lubricated packings. 

O 



INTRODUCTION 

-^he mechanical properties of granular packings are 
L_^nsitive to fine geometric details of the particle ar- 
rangement, as very small motions can modify the 
^^orce-carrying contact network. In practice, direct 
i/Tneasurements of such important internal state vari- 
O^bles as the coordination number and the distribution 



ers ([Lings & Greening 2001 ). The recent soil 



mechanics literature ( [Thomann & Hryciw 1990 



^^f contact orientations ( [Radjai & Roux 2004[ ) are usu- 
ally impossible. 

O In this context, ultrasonic wave propagation mea- 
^^urements might provide useful information on the 
Vintemal structure of granular samples under confin- 
.^ng stresses. Wave propagation has been used as a 
^<^robe to investigate the micro structure of granular 
packings by several physics groups (I Jia et al. 19991 
■ - %a & Mills 200n IGilles & Coste 20031) . At low fre- 
quencies such that the wavelengths are very long 
compared to the heterogeneity of the medium, the 
granular medium is effectively a homogeneous con- 
tinuum to the propagating wave, while at high fre- 
quencies when the wavelength decreases down to the 
order of the grain size, scattering effects caused by 
the spatial fluctuations of force chains lead to diffu- 
sive transport of sound waves (jJia 2004|) . 

Wave propagation has also become a stan- 
dard method to measure elastic moduli in 
geotechnics laboratories ( [Geoffroy et al. 2003 ), 
where rheological testing devices are often 
equipped with specially designed transduc- 



Geoffroy et al. 2003) made it clear that "dynamic" 
measurements of elastic moduli (wave propagation 
or resonance modes) agree with "static" ones (slopes 
of stress-strain curves), provided strain increments 
are small enough (below 10^^). Correctly measured 
elastic moduli therefore determine long wavelength 
sound velocities in granular materials as in ordinary 
solids. 

Discrete particle simulations can also be used to 
evaluate elastic properties of model granular mate- 
rials (IRoux 19971 IMakse et al. 20041 IRoux 2004 ). If 
experimental packings are correctly simulated, such 
studies can clarify the relations between wave propa- 
gation measurements and sample microstructure. 

We first briefly report here (Sec. O on experimen- 
tal measurements on pressure-dependent sound veloc- 
ities and attenuation in dense samples of dry glass 
beads (hereafter denoted as El samples), as well as 
lubricated ones in which the effects of intergranular 
friction are strongly reduced in the preparation stage. 
The resulting materials (denoted as E2 samples) are, 
remarkably, less dense, but stiff er (with notably larger 
elastic wave velocities). 

In order to investigate the microscopic origins of 
these results, we use numerical simulations, as pre- 
sented in Sec. [3l which is the main part of the present 
communication. We resort to suitable, yet simplified. 
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models to mimic laboratory assembling procedures 
(thus producing sets of samples denoted as A and B 
by a lubricated procedure, and a third series C by a 
"vibrated" procedure, see Sec. B.ll for details). Com- 
parisons of numerical results to experimental ones 
lead to an interpretation of the observed differences 
between dry and lubricated samples in terms of coor- 
dination number. 

2 EXPERIMENTS 

The glass beads used in our experiments are of diam- 
eter d = 0.3-0.4 mm, randomly deposited by pouring 
and vibrating in a duralumin cylinder of diameter W 
= 30 mm and varying height from 10 mm to 30 mm. 
The container is closed with two fitting pistons and a 
normal load is applied to the granular sample across 
the top piston. Before the ultrasonic measurements, 
one cycle of loading-unloading is performed in the 
granular packing in order to consolidate the sample 
and minimize its hysteretic behavior. A plane-wave 
generating transducer of diameter 30 mm (top pis- 
ton) and a detecting transducer of diameter 30 mm 
are placed on the axis at the top and bottom of the 
cylindrical container in direct contact with the glass 
beads. Once the cell is filled, a normal load corre- 
sponding to apparent pressures P ranging from 30 kPa 
to 1000 kPa is applied to the upper piston using a 
jackscrew arrangement, while the lower piston is held 
fixed, as in oedometric loading. (The preloading pres- 
sure cycle reached up to 400 kPa). Broadband short 
pulse excitations of 2 jis duration centered at a fre- 
quency of 500 kHz are applied to the source trans- 
ducer. The time of flight of the transmitted ultrasonic 
signal is used to measure the sound speeds, which are 
shown on figure [H and to which results of numerical 
simulations will be compared in Sec. 13.21 (Figure |3]). 
Measurements are first performed on assemblies of 
dry beads, initially prepared in a dense state by lay- 
erwise deposition and tapping (samples El). Then, to 
study the influence of intergranular friction on the ini- 
tial structure and elastic properties of the medium, we 
mix a small amount of liquid lubricant (trioleine, vol- 
ume fraction 0.5%) with the granular sample for tens 
of minutes to distribute the oil uniformly among the 
grains. We observe (figure [U), that the solid fraction 
in the obtained lubricated samples (denoted as E2) 
is about 0.62, lower than in the dense dry ones El 
( ~ 0.64), while sound wave speeds are significantly 
higher for E2 specimens (by 10 to 20%) than in the 
El case. It is also apparent on Figured] (a logarith- 
mic plot), that the increase of sound velocities with 
pressure is slightly faster in El samples. 

3 NUMERICAL SIMULATIONS 

We now turn to numerical simulations to explore the 
origins, at the scale of the contact network, of the dif- 
ferences between the dry (El) and the lubricated (E2) 
samples. 
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Figure 1 . Velocity of longitudinal (bottom plot) and transverse 
(top plot) sound waves as a function of confining pressure for 
laboratory samples El (dots, assembled by vibrating dry grains) 
and E2 (crosses, lubricated). 

3.1 Sample preparation and characterization. 

Our numerical simulations, like oth- 
ers (Maks e et al. 20041) (but, admittedly, unlike 
the experiments of Sec . II) focus on homoge- 
neous, isotropic states, under varying pressure 
P. The samples - the same as those studied 
by 'Agno lin & Roux (2005 [ )- comprise 4000 identical 
beads of diameter a in a periodic cell which changes 
size as stresses are applied. Averages and standard 
deviations are evaluated on 5 such samples. 

Spherical grains are attributed the elastic proper- 
ties of glass beads (Young modulus = 70 GPa, 
Poisson coefficient v = 0.3), and assembled by a 
standard molecular dynamics method. The contact 
laws involve Coulomb friction and a suitably sim- 
plified form of (Hertz-Mindlin) contact elasticity 



(see ( Agnolin & Roux 2005 1 and references therein). 

The most frequently used numerical proce- 
dure (IMakse et al. 20 04) to obtain dense configura- 
tions is to suppress friction altogether while com- 
pressing a granular gas to mechanical equilibrium un- 
der a given isotropic pressure, thus simulating per- 
fect lubrication. This results in typical isotropic ran- 
dom close packing structures (lO'Hern et al. 2003)) . 
the properties of which are not sensitive to the de- 
tails of the procedure, provided it is fast enough to 
bypass crystal nucleation entirely. On applying this 
method with prescribed pressure P = 10 kPa, we ob- 
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tained samples we denote as A. Their coordination 
number z*, evaluated on eliminating inactive grains 
("rattlers") from the count, is close to 6, the value 
it should approach (IRoux 20001 IMakse et al. 20041) in 
the rigid limit, as the average contact deflection h be- 
comes negligible, h/a(x {P/ E^/^ 0. On imposing 
larger pressure levels (up to quite high values in sim- 
ulations) a friction coefficient fi = 0.3 is introduced. 
This implicitly assumes that the perfect lubrication of 
the assembling stage (in which intergranular forces 
are transmitted through a very thin layer of lubricant) 
disappears as large static pressures bring solid sur- 
faces into contact. Another series of numerical sam- 
ples, type B ones, are also made assuming only a very 
low friction in the initial stage (/io = 0.02), as a first 
model of imperfect lubrication. Finally, in order to 
imitate the preparation of dense dry samples of type 
El in the laboratory, a third series C of initial states 
under 10 kPa is prepared as follows. First, A config- 
urations are slightly dilated, scaling coordinates by a 
common factor A = 1.005; then grains are mixed at 
constant volume, as by thermal agitation : the sys- 
tem is thus strongly shaken in a dense state; the fi- 
nal step is a compression, with friction (/i = 0.3) and 
some viscous dissipation to mechanical equilibrium 
at P = 10 kPa. 

Figure [2] displays the values of solid volume frac- 
tion $ and coordination number z* as functions of 
pressure P in states A, B and C. While B samples are 
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Figure 2. Volume fractions and coordination numbers as func- 
tions of P in isotropically compressed samples of types A 
(square dots), B (stars), and C (open circles). 

similar to A ones, with somewhat lower values of $ 
and z*, C configurations, remarkably, although very 
nearly as dense as A ones (as expected given their 
preparation method), and actually denser than B ones, 
exhibit much lower coordination numbers. The differ- 
ence between C and A states is gradually reduced as 
P grows. 

3.2 Wave velocities. 

Isotropic states A, B, C possess two independent elas- 
tic constants, the bulk (B) and shear (G) moduli, from 
which the velocities of longitudinal (Vp) and trans- 



as Vp 



'B + 4G/3 



and Vs = \ l —, P denoting the 

mass density in the granular material. 

Once mechanical equilibrium states are obtained 
with sufficient accuracy, the stiffness matrix of the 
contact network is built, and both elastic moduli are 
obtained on solving a linear system of equations, in 
which particle displacements and rotations as well as 
strain increments are the unknown, and the imposed 
stress increments determine the right-hand side. Fig- 
ure [3] presents the simulated values of Vp and for 
all three sample types and compares them to exper- 
imental results. Scarcely coordinated C samples are 
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Figure 3. Comparisons between numerical and experimental 
values for Vp (bottom plot) and Vs (top part). Numerical re- 
sults are plotted as connected symbols, as on figure |2l square 
dots for A (perfect lubrication), stars for B (imperfect lubrica- 
tion), open circles for C ("vibrated" dry grains). Experimental 
results are shown as filled circles (El, dry samples) and crosses 
(E2, lubricated samples) for P in the 70 kPa-1 MPa range, as 
on Figure [T] while the data obtained by Domenico on dry glass 
beads are plotted as black dots for P around 10 MPa. 

clearly in much better agreement with experimental 
data on dense, dry packings than A or B ones in the 
P ~ 100 kPa range. The difference between B and 
C systems is qualitatively the same as the one be- 
tween dry and lubricated laboratory packings: a lower 
density, but higher wave velocities that increase a lit- 
tle slower with pressure. Some indications about the 
pressure depence of elastic moduli B and G in A 
and C configurations, and a discussion of their pre- 
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diction by micromechanical modeling schemes, are 



given by Agnolin & Roux (2005 1 in these proceed 



ings. Despite the crudeness of the numerical models 
for "shaking" (C) or "lubrication" (A, B), they appear 
to capture the right experimental trends. We suggest 
therefore that the larger stiffness (or sound speed) in 
lubricated packings is due to their larger coordination 
number. 

We therefore conclude that C-type samples are bet- 
ter models for dense dry granular packings (of type 
El). Such a statement apparently contradicts the good 
agreement reported by Mak se et al. (2004] ) between 
numerical measurements of sound speeds on sam- 
ples of type A and experiments on dry specimen of 
type El, such as those by [Domenico (1977| ), shown on 
Figure [3l However, on confronting simulation results 
with experimental ones, those authors focussed on a 
higher pressure range (above several MPa, see Fig- 
ure [3]), in which A and C-type systems are much less 
differentiated. We checked that our numerical results 
on A samples coincide nearly perfectly with those of 
Makse et al. (2004 ]). Thanks to our experimental re- 
sults for smaller confining pressures we can also dis- 
tinguish between dense systems with high and low co- 
ordination numbers, the latter ones being more appro- 
priate as models for dense, dry granular assemblies. 

It might also be noted that in a different communi- 
cation in the present proceedings (' Roux 20051) . it is ar- 
gued that the mechanical properties of C-type samples 
in quasistatic, axisymmetric triaxial compression are 
closer to those observed in the laboratory with sand or 
glass bead samples assembled by pouring and shaking 
than those of A configurations. 

4 CONCLUSION 

Numerical simulations show that wave propagation 
in dense bead packings assembled with various pro- 
cedures can reveal differences in their microstruc- 
ture, beyond the sole density. Admittedly, more ac- 
curate comparisons with experiments are certainly 
necessary. This requires a more detailed knowledge 
of laboratory samples and their anisotropic state of 
stress (|Khidas & Jia 20051) . a possible consideration 
of the effects of capillary adhesion, as well as more 
realistic numerical models for laboratory procedures 
(|Emam et al. 20051) . However, the results presented 
here show that discrete numerical simulations can re- 
late experimentally accessible data on wave speeds 
to internal variables such as coordination number or 
fabric. They clearly stress the need for a better un- 
derstanding of the influence of the preparation pro- 
cedure on the subsequent mechanical properties in 
quasistatic conditions. They suggest a plausible in- 
terpretation of the observed larger wave velocity in 
lubricated systems. Numerical procedures, even for 
isotropic samples, should not be selected as appro- 
priate because they produce the right density. Coor- 
dination number, a largely independent state variable. 



strongly influences elastic properties and quasistatic 
stress-strain laws (|Roux 2005 1) . 
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